# from geopy.distance import geodesic
# 导入核密度估计包
import pandas as pd
from pykrige.ok import OrdinaryKriging
import numpy as np
from math import radians, cos, sin, asin, sqrt
# import numba
from numba import jit, njit


@njit
def haversine(lon1, lat1, lon2, lat2):
    R = 6372.8
    dLon = radians(lon2 - lon1)
    dLat = radians(lat2 - lat1)
    lat1 = radians(lat1)
    lat2 = radians(lat2)
    a = sin(dLat / 2) ** 2 + cos(lat1) * cos(lat2) * sin(dLon / 2) ** 2
    c = 2 * asin(sqrt(a))
    d = R * c
    return d


def IDW(x, y, z, xi, yi):
    lstxyzi = []
    for p in range(len(xi)):
        # lstdist = []
        nearlist = []
        orign = None
        for s in range(len(x)):
            d = (haversine(x[s], y[s], xi[p], yi[p]))
            # d = geodesic(kp, latlon).m
            if d == 0:
                orign = z[s]
                break
            elif d < 10:
                # lstdist.append(d)
                nearlist.append((d, z[s]))
        if orign is not None:
            u = orign
        else:
            # nearlist.sort()
            # nearestArr = np.array(nearlist[:30])
            if not len(nearlist):
                u = np.nan
            else:
                nearestArr = np.array(nearlist)

                sumsup = list((1 / np.power(nearestArr[:, 0], 2)))
                suminf = np.sum(sumsup)
                sumsup = np.sum(np.array(sumsup) * nearestArr[:, 1])
                u = (sumsup / suminf).round(2)
        xyzi = [xi[p], yi[p], u]
        lstxyzi.append(xyzi)

    return (lstxyzi)

def krige(lons, lats, data, grid_lon, grid_lat):
    # OK = OrdinaryKriging(lons, lats, data, variogram_model='gaussian')
    OK = OrdinaryKriging(lons, lats, data, variogram_model='linear')
    # start_lat, end_lat, lat_reso, xdim = 28.0, 19.0, -0.01, 900
    # start_lon, end_lon, lon_reso, ydim = 102.0, 114.0, 0.01, 1200
    # grid_lat = np.linspace(start_lat, end_lat - lat_reso, xdim)
    # grid_lon = np.linspace(start_lon, end_lon - lon_reso, ydim)  # For shape compatibility
    z1, ss1 = OK.execute('grid', grid_lon, grid_lat, backend='loop')
    # z1, ss1 = OK.execute('grid', grid_lon[:920], grid_lat[:690])
    # print("===>", z1.data.shape, z1.max())
    # 将插值网格数据整理
    # xgrid, ygrid = np.meshgrid(grid_lat, grid_lon)
    # df_grid = pd.DataFrame(dict(long=xgrid.flatten(), lat=ygrid.flatten()))

    # 添加插值结果

    # df_grid["Krig_gaussian"] = z1.flatten()
    return z1


if __name__ == '__main__':
    haversine(102.01, 25.01, 102.02, 25.02)
